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ABSTRACT 

A comparison of standard z transform and algebraic 
substitution synthesis methods for digital filters are 
presented with emphasis on the frequency domain character- 
istics. A generalization of Martin's procedure, which leads 
to recursive filters, is obtained. Also a new method to 
synthesize a digital filter from a magnitude squared versus 
frequency characteristic specification is presented. 
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I. 



INTRODUCTION 



Digital filters have become increasingly important 
because of technological developments in integrated circuits 
and digital techniques. This progress has increased the 
number of applications in which digital filtering techniques 
can be used. In spite of the fact that it is a relatively 
new development, a vast body of literature is available for 
digital filters. The relationship between digital filters 
and classical continuous filter theory can be depicted as 
in Fig. 1.1. 

Start from the differential equation at the center of 
the diagram. The Laplace transform of the differential 
equation may be taken and a transfer function formed. The 
frequency spectrum of the continuous filter may be evaluated 
as shown in the diagram by letting s = jco. The upper part 
of the diagram shows the relationship between the transfer 
function of the continuous filter in the frequency domain, 
and the frequency domain characteristics of the sampled impulse 
response of the continuous filter. In deriving these results 
use is made of the Z transform. The lower part of the diagram 
shows numerical integration of the differential equations. 
Adopting a numerical integration algorithm yields a difference 
equation which can generally be expressed as a discrete time 
equation. For a given integration scheme, as illustrated in 
the diagram, it is possible to substitute a function of z 
(depending upon the type of integration) for the variable s 
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Fig, 1,1. Relationship between digital filter and continuous 
filter theory. 



of the Laplace transformed continuous equation. When z = 
exp (jwt) is substituted in this equation a spectrum is obtained 
which, in general, is not identical to the spectrum of the 
sampled impulse response of the continuous filter. One of the 
objectives of this thesis is to compare numerical integration 
with the standard Z transform results as depicted in the upper 
portion of the diagram. 

The Z transform technique is described in the work of J. 
Ragazzini and L. H. Zadah [1] and has been discussed by many 
other authors including B. C. Kuo [2]. J. Salzer [3] has 
written a classic paper on the comparison of digital integra- 
tors in the frequency domain with the ideal continuous 

2 z-1 

integrator. The bilinear substitution, s = — j - , which 
corresponds to trapezoidal integration, was considered very 
early by A. Tustin [4] and has been discussed thoroughly by 
K. Steiglitz [5] . 

In this thesis several areas of digital filter theory are 
considered in detail. 

Chapter two contains a derivation of basic Z transform 
theory, a review of the synthesis of digital filters using 
standard Z transforms, and a discussion of the effect of 
sampling frequency on the frequency spectrum. The relationship 
between the s and z plane representation of the equations is 
presented. Chapter three contains a general interpretation of 
the algebraic substitution, s = f(z), for several numerical 
integration techniques, a discussion of the bilinear substitu- 
tion (trapezoidal integration) , and a development of the effects 
of sampling time on the precision of numerical calculations. 
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Chapter four discusses the synthesis of digital filters from 
frequency spectrum characteristics. It contains a derivation, 
discussion, and a practical application of a direct method 
of synthesis originated by Martin [6] . A new generalization 
of this method to obtain recursive filters from a specification 
of magnitude versus frequency characteristics, is presented. 

The advantages of this technique include the fact that a zero 
phase shift filter can be obtained. Derivation of another new 
method of finding the coefficients of a digital filter, starting 
from a magnitude squared criterion, is presented. The advan- 
tages of this technique include reducing computer storage 
requirements, and reducing the number of multiplications to one 
half of those required by the Martin method. Chapter five 
presents a summary and several questions for future research. 
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II. THE Z TRANSFORM TECHNIQUE APPLIED 
TO DISCRETE TIME SIGNALS AND DIGITAL FILTERS 



Systems may be classified as linear, time invariant, 
causal and so on. They may also be classified according to 
the nature of the signals present. A continuous signal is 
expressed as a continuous (or piecewise continuous) function 
of the independent variable, t. The signal must be uniquely 
defined at all values of t within a given range, except 
possibly at a denumerable set of points, as for example in a. 
square wave. A discrete signal is defined only at a sequence 
of discrete values of the independent variable t. A quantized 
signal is one which can assume only a denumerable number of 
different values, but quantized signals may be either continuous 
or discrete, as discussed by P. M. DeRusso et al , [7]. 

There are many examples of continuous systems in nature, 
and they are to be found almost everywhere. Examples of 
discrete and quantized systems are perhaps more difficult to 
conceive but this does not mean they are scarce.. For instance, 
the output of a continuous system when sampled for digital 
computation, or other purposes, becomes a discrete system. 

Some systems are discrete by nature. The distance to an object 
as measured with a pulsed radar; or a pulse coded communica- 
tions system which is time multiplexed are examples of inherently 
discrete physical systems. A digital computer is an excellent 
example of a discrete, quantized system. The quantizing levels 
are determined ultimately by the finite number of bits of the 
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number representation. The discrete characteristic is deter- 
mined by the fact that the independent variable assumes values 
at finite increments only. 

In order to deal with discrete time quantized systems, 
existing Z transform theory has been extended [Ref. 2] .. This 
theory was developed primarily for sampled data control 
systems and resembles Laplace transform theory for continuous 
systems in that it is operational in nature. The Z transform 
has proven to be useful for the study of all discrete time 
systems and it is extended approach which is reviewed in this 
chapter. Several important fundamentals are developed and a 
basic misapplication of the theory which occurs commonly in 
the literature, is discussed. 

V 

A. DEFINITION OF THE Z TRANSFORM OF A DISCRETE SIGNAL 

Consider a continuous signal x(t), for t>0, which is 
sampled at a frequency, f or equivalently, every T = 1/f 

S o 

seconds . 

Assume that the signal is sampled ideally, that i s, the 
sampled signal assumes the value of the continuous signal with 
infinite accuracy at discrete values of the independent 
variable, t=nT;n= 0,1,2, ... 

The process of sampling can be described as 

x* (t) = x ( t ) s ( t ) (2.1) 

where 

x* (t) = the sampled signal 
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x(t) ~ the continuous signal 

s(t) = the sampling function, representing a train of 



impulses, T seconds apart. 



s(t) = Z 6 (t-nT) 
n=0 



From (2.1) and (2.2) 



x* (t) = Z x(nT) 6 (t-nT) 
n=0 



Taking the Laplace transform of the sampled signal 
have from Eq . (2.3) 

X* (s) = £ (x* (t) } 

oo 

= Z x(nT) £ (6 (t-nT) } 
n=0 

But £ (6 (t-nT)} = e _nTs 
so that Eq. (2.4) becomes 

OO 

X* (s) = Z x(nT) e“ nTs 
n=0 

Now, define a new variable z, such that 
sT 

z = e 

so that Eq. (2.6) becomes 



X (z) = X* (s) 



z=e 



= Z x(nT) z 
sT n=0 



-n 



The standard procedure to find the Z transform of 
discrete signal can be summarized as follows: 



( 2 . 2 ) 

( 2 . 3 ) 

we 

( 2 . 4 ) 

( 2 . 5 ) 

( 2 . 6 ) 

(2.7) 

( 2 . 8 ) 
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(a) Sample the continuous signal 

(b) Find the Laplace transform of the sampled function 

sT 

(c) Substitute e = z in the Laplace transform of the 
sampled signal to obtain the Z transform. 

Using common notation, the above can be written as follows: 



X (z) = ^ {x (t) } 



[£{x* (t) }] 



z=e 



sT 



(2.9) 



This process of obtaining the Z transform of a signal is 
illustrated in Fig. 2.1. 

1 . Obtaining the Z Transform in Closed Form 

The foregoing definition of the Z transform provides 
a method for obtaining the transform of a signal, but it has 
two disadvantages: 

a. Results are presented in open form because the 
transform of a signal is given as a summation of inverse 
powers of z whose coefficients are the values of the signal at 
the sampling times. For sampled functions derived from contin- 
uous functions with poles in the left half s plane, such 
coefficients will eventually decay to zero. If the decaying 
process is slow compared to the sampling interval, a large 
number of terms in the sampled Z transform expansion are required 
to represent the successive values. 

b. The process of obtaining the Z transform may 
be very long, particularly when only the Laplace transform of 
the continuous signal is known. The process requires taking 
the inverse Laplace transform of the continuous signal, sampling 
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Procedyre to obtain the z Transform 




the resulting time function, and finally expressing the sample 
values as a series of inverse powers of z. 

There is an alternate method for obtaining the Z 

transform of a signal that is particularly suitable when the 

Laplace transform of the continuous signal is known. This 
latter approach gives the result in closed form. From Eq. 
( 2 . 1 ) 

X* (s) = £ (x(t) s (t) } (2.10) 

and so X* (s) = x (s) * S(s) (2.11) 



where the asterisk denotes convolution. 

c+j«> 

X* (s) = 2 ~j j X (X) S(s-X) dX 

C- j°° 



aOO CO 

now S(s) = E 6 (t-nT) e st dt 
J o n=0 



S (s) 



= E 
n=0 



-nTs 

e 



( 2 . 12 ) 



S(s) = 1 / (1-exp (-s T) ) (2.13) 

provided |exp(-nT)|< 1 

In (2.11), use has been made of the sifting integral 
property of the impulse, namely 

6 (t-nT) f(t) dt = f (nT) (2.14) 

■* o 
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Therefore, the Laplace transform of the sampled 
function can be written in terms of the complex convolution 
integral 



X* (s) 



1 

2tt j 



c+ j 00 
XU) 
c- j°° 



1 



1-e 



- (s-X) T 



dX 



(2.15) 



If as shown in Fig. 2.2 is taken as the path 
enclosing the left half of the s plane then integral (2.15) 
may be divided into two integrals as follows 



X* (s) 



— 4 *(*) dX 

2 7T j •'r i 1- e 1 ’ 



j 

2n j ' 



X(X) 



1- e 



- (s-X) 



dX 



(2.16) 



Note that is composed of the straight line joining 
the points (c,-j°°) and (c,+j°°) and the counter-clockwise semi- 
circle or infinite radius. 

If the condition 

lim s X(s) = 0 (2.17) 

s 00 

is met, then the second integral tends towards zero and Eq. 
(2.16) becomes 

X* (s ) = x(X) 1 , , dX (2 .18) 

2ir j ^ 1- e" (S_A) 

A very important case occurs when the Laplace 
transform X(s) of the continuous signal is a ratio of two 
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The contour integration path 
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two polynomials 



X(s) 



N(s) 

D(s) 



(2.19) 



In this case the function, X(s), has a finite number 
of poles and the integral around the closed path can be 
evaluated using the residue theorem. Thus Eq . (2.18) becomes 



X* (s) = E 



N(X n ) 



n=l D' (X ) 1- e 
n 



-(s-A n )T 



( 2 . 20 ) 



where k is the number of simple poles of X(A) located at A^, 
with 



, ( X ) = d D( s ) 

n ds 



( 2 . 21 ) 



3 = X 



n 



and substituting Eq. (2.7) into (2.20), yields X(z) the Z 
transform of x(t) 



k N ( A ) , 

X (z) = Z t — m ‘ (2.22) 

n-1 D ' (X ) 1- e n z 
n 

The above formula (2.22) represents a closed form 
for the Z transform of the signal which contains as many terms 
as the degree of the denominator polynomial of the function 
X(s). Expansion of Eq. (2.22) in a series of inverse powers 
of z yields Eq. (2.8). 

B. SYNTHESIZING A DIGITAL FILTER WITH THE STANDARD Z TRANSFORM 
In the previous section two methods have been developed 
to obtain the Z transform of a sampled signal, starting from 
either the continuous time function or its Laplace transform. 
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The same procedure can be followed to obtain the Z transform 
corresponding to the sampled output of a continuous filter, 
considering the signal to be the impulse response of the 
filter, h (t) . 

In general there are many known procedures for the design 
of continuous filters to meet varioas types of specifications. 
These procedures are usually described in the s domain or in 
the frequency domain, s = jto. The Z transform is used to 
translate continuous time realizations into discrete time 
realizations or digital filters. 

In control theory, when the differential equations of 
the systems are known, it is often necessary to simulate 
continuous systems on a digital computer. One approach is to 
consider the impulse response which is the inverse Laplace 
transform of the system's transfer function.. The Z transform 
can be applied to this signal after sampling, yielding a discrete 
simulation of the system. 

Equation (2.8) assures that the Z transform of a digital 
filter or system, when expanded in terms of inverse powers 
of z, represents identically the impulse response of the 
continuous filter or system sampled at intervals T. 

If X( z) is the Z transform of the input signal to a system, 
and if H(z) is the Z transform of the system transfer junction, 
then 

Y (z) = X (z) H (z) (2.23) 
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Here Y (z) represents the Z transform of the output signal. 



and corresponds to the numerical convolution of the input 
signal an the system transfer function, as is demonstrated 



below 




X(z) 


oo 

= Z x(mT) Z -IT1 

m=0 


H ( z ) 


co (2.24) 

- E h(kT) z k 

k=0 



where h(kT) is the sampled impulse response of the system. 
Substituting (2.24) into (2.23) 



Y(z) 


00 OO 

= E E x(irtT) h (kT) z~ (k+m) (2.25) 

m=0 k=0 



and replacing subscripts, such that k+m - n or k = n-m, then 
the limits of the summations are converted as follows: 
k = 0 becomes n = m and k = 00 becomes n = °°.. 



Thus 


(2.25) becomes 


Y(z) 


OO oo 

= Z Z x(mT) h ( (n-m) T) z” n (2.26) 

m=0 n=m 



Changing the order of the summations and neglecting or 
adding zero terms of the type h(-kT) (k=l,2, ) yields 



Y(z) 


°o n 

= Z Z X(mT) h ( (n-m) T) z n (2.27) 

n=0 m=0 




OO 

= Z y(nT) z n 

n=0 
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Hence the coefficient y (nT) of z n in Y(Z) is given by 
n 

y (nT) = E x (nT) h((n-m)T) (2.28) 

m=0 

Equation (2.28) is the numerical convolution of the input 
samples and the digital filter sampled impulse response. 

The output signal transform, Y(z), when expanded in a 
series of inverse powers of z represents the output signal at 
different times. If X(z) and H(z) are obtained from sampling 
a continuous input signal and the impulse response of a con- 
tinuous system or filter respectively, then the coefficients 
in the expansion of Y(z) will correspond to the sampled output 
of the continuous system or filter. Hence when synthesizing 
a digital filter using the Z transform of a continuous filter 
impulse response sampled at intervals T, the resulting output 
signal is an exact description in the time domain in the sense 
that the output values of the digital realization will be equal 
to the values of the continuous filter output signal sampled 
at intervals T. 

Regarding the frequency spectrum of the filters, consider 
a continuous system with transfer function H(s) . If an input 
signal E^(s) is applied, the output signal E q (s) has a Laplace 
transform 

E q (s) = H (s) E i (s) (2.29) 

The inverse Laplace transform of E q (s) is the output 
signal e Q (t) = £ ^{E o (s)}, which now is sampled at the rate of 
1 / T samples per second, yielding the sampled output signal 
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e Q * (t) as shown in Fig. 2.3a. The Fourier transform of this 
signal is the desired output spectrum 

E * (ju) = 3 {e *(t)} (2., 30) 

o o 

This spectrum could also have been obtained by convolving 

the Fourier transform of the sampling function ? 6 [j (to — 

n=-°° 

nto ) ] such that 
s 

oo 

E *(jto) = E ( j to ) * E 6 [ j (to-nto )] (2..3I) 

° ° n =-» . S 

CO 

= £ E o ( j [tt-nco ] ) 

n=-°° 

This signal processing may be compared with the develop- 
ment of the digital filter shown in Fig. 2.3b. We assume the 

same sampling frequency f = 1 / T and obtain the Z transform 

s 

H(z) of the filter function H(s). Then the Z transform of the 
filter output E q (z) is given by 

E (z) = H (z) E ± (z) ' (2.32) 

where E^(z) is the Z transform of the sampled filter input 
e^* (t) . Substituting z = exp(jtoT) into E 0 ( z ) yields the desired 
output spectrum E Q *(jto). The following question may now be 
asked: Is this output spectrum of the digital filter identical 

with the output spectrum obtained in Eq. (2.31)? This question 
is discussed in the next section. 

1 . The Effect of Sampling on the Frequency Spectrum 

For a signal of finite bandwidth w radians called a 
band limited signal shown in Fig. 2.4a Shannon's Sampling 
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Fig. 2.3. Equivalence between discrete and continuous 
filters at sampling instants. 



Theorem [8], requires a sampling rate equal to at least twice 
the maximum frequency w of the signal in order to recover the 
signal without distortition after passing it through an. ideal 
low pass filter. The spectrum of the signal sampled at a 
higher rate than the Nyquist rate is shown in Fig. 2.4b.. If 
the sampling frequency is lower than the limit set by Nyquist, 
then some overlapping of spectrum will occur. This is known 
as the aliasing effect and is illustrated in lig. 2.4c.. It has 
been shown by L. J. Fogel [9], that, it is possible to recover 
a signal sampled at a lower rate than Nyquist (or one for which 
aliasing has taken place) if for every sample additional infor- 
mation about the signal is known. For example, this is possible 
if its derivatives with respect to time are known. However, 
this idea cannot be applied here because we restrict ourselves 
to single input filters where only the values of the signal are 
known . 

To answer the previous question regarding comparison 
of the two spectra of Fig. 2.3, consider the continuous case 
of Fig. 2.3a, where for s = jw 

E Q (jo)) = H(jw) E i (jto) (2.33) 

After sampling this spectrum, E q * ( jw) is the convolu- 
tion of E Q (jw) and the spectrum of the sampling signal, which 
is a train of impulses separated by w s on the w axis as given 
by Eq . (2.31) . The effect is to make the sampled spectrum 

repetitive with period u) g . From Shannon's Sampling Theorem, 
the sampled output can be reconstructed without aliasing distor- 
tion provided E Q (ja)) is band limited and the sampling frequency 
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(a) Spectrum of a band limited signal 



|E ( jw) | 

w 2 ti f 

s 



(b) Spectrum of a sampled band limited 
signal, f > 2w 



|E ( jw) | 

l 

— V/ — ! vr 



L ^-» 2TTf 0) 

s 

aliasing 

effect 



(c) Spectrum of a sampled band limited 
signal, 2irf s < 2w 



Fig. 2.4. Effect of sampling frequency on the 
frequency spectrum of the sampled 
signal . 
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is higher than two times the maximum frequency of E o (jw). For 
E (jw) to be band limited either E^(jm) or H(jw) have to be 
band limited. 

Consider now the discrete case shown in Fig.. 2.3b, 
where e^ (t) is sampled and the Z transform is obtained. From 
Eq. (2.8) 



E i (z) 



E i * (s) 



z 



sT 



e 



(2.34) 



V (e 



Substituting Eq. 
^“ T ) = H* (e jwT ) 



(2.34) into equation (2.32) yields 
E i *(e jwT ) (2.35) 



1 tgjT 

As may be seen, the output spectrum E q * (e J ) is the 

product of the two spectra H* (e-* wT ) and E^* (e^ uT ) .. Clearly 

both these spectra must be undistorted (no aliasing) if 

E * (e- |wT ) is to be undistorted. Hence for this case the 
o 

requirement is different from the case of Fig. 2.3a.. In view of 
the above remarks we conclude that it is necessary to. have 
E^(ju) and H(jm) band limited and furthermore, one must sample 
at a frequency above the Nyquist rate for both spectra of 
Figs. 2.3a and 2.3b to be identical. 

As an example consider the design of a digital low 
pass filter with a cut off frequency of 2 KHz which is required 
to filter a signal of 10 kHz bandwidth. Such a filter must 
be designed with a sampling rate greater than 20 kHz.. Note 
that if the signal had been obtained from a continuous filter, 
it could be recovered by sampling at a rate of only 4 kHz. 

It is clear from the above remarks and the example given, 
that there is a maximum frequency that can be handled' by a 



25 



digital filter in real time operation. In fact this limit will 
depend on the computing speed of the processor and the general 

complexity of the filter. It is desirable to keep digital 

filter designs as simple as possible, although considerations 
of stability may modify these considerations. 

2 . The Interelationship Between the s and z Planes 

The interelationship between the s plane and the z 

plane can be derived from the definition of the variable z 
namely 

z = e sT (2.36) 

A point in the s plane, namely (a+jio) will be 
mapped into the point in the z plane 

2 = e (o+ jw) T _ e 0T e j M T (2.37) 

and so I z | = e aT 

Therefore, the points in the s plane with constant 

damping, i.e., a = constant, map onto a circle of constant 
0 T 

radius, e , in the z plane. Points in the s plane with 
constant frequency, i.e. w = constant, map onto radial lines, 
arg z = toT. 

Furthermore, for points in the left half s plane 
o < 0 , and so 

|z| = e aT < 1 for T > 0 (2.38) 

Equation (2.38) indicates that the entire L-E..P. 
maps into the inside of the unit circle in the z plane. 
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Although useful, the foregoing discussion may be 



misleading unless it is applied properly. The misuse of 
Eg. (2.36) occurs often in the literature. Take for example 
the paper by C. J. Greaves and J. A. Cadzow [10] . Their mis- 
take consists of replacing the variable s in the Laplace 
transform of the continuous filter by its corresponding z, with 
s = (1/T) In z as 



H ( z) = H (s) 



(2 ..39a) 



s = (1/T) In z 



This equation is in error and should be replaced by 



H (z) 



H* (s) 



s 



(1/T) In z 



(2. .3 9b) 



The reason is clear from the closed form expression 
for the Z transform, Eq. (2.18), where before the substitution 
of variables is made, a complex integral has to be evaluated. 

However, there is a direct connection between the 
Laplace transform of the continuous function and the Z trans- 
form of the sampled function. Namely the poles of the Laplace 
transform of the continuous function are mapped into the z plane 
of the sampled function using the mapping function of (2.36). 

This can be seen from equation (2.22) where every 
term in the summation will have a pole in the z plane at 

z = e^ nT (2.39c) 



which indicates the mapping of the singularity A n according to 
the mapping relation of (2.36). 
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In summary Eq . (2.36) may only be used to transform 

poles of the continuous function in the s plane to pole 
locations of the sampled function in the z plane. It does not 
represent a change in variable except when applied directly 
to a sampled signal. As noted this subtlety has proven to be 
a source of many mistatements in the literature. 
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III. THE ALGEBRAIC SUBSTITUTION METHOD 



An equivalent digital representation of a continuous 

system or filter is often required in practice. This could 

arise, for example, in a digital simulation of a continuous 

control problem, or in an attempt to extend continuous filter 

theory to provide digital filter realizations. This latter 

application is very appealing because continuous filter theory 

is very well developed and procedures exist to match given 

specifications to filter relizations. Designs have been 

tabulated and are available in many handbooks. If a direct 

method to obtain the equivalent digital realization of a contin- 

/ 

uous filter exists, it would make all previous knowledge of 
continuous filter theory available for the discrete case.. 

It is convenient when implementing a digital filter, to 
obtain its representation as a Z transform, and particularly 
as a ratio of two polynomials of inverse powers of z. Such a 
representation indicates the linear combinations which have to 
be performed upon past values of input and output signals 
samples, to obtain the present value of the output. 

The most common description of continuous filters in the s 
domain is as a ratio of two polynomials in s and so it is con- 
venient to have a rational function of z to directly replace 
for s in the continuous type description. This is if 

s = f(z) ((3.1) 

then H (s) = H(f(z)) = G(z) . (3.2) 
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If the function f(z) is a ratio of two polynomials of 
inverse powers of z and if H(s) is a ratio of two polynomials 
in s, then G(z) in Eq. (3.2) is assured to be the ratio of two 
polynomials in inverse powers of z. 

This procedure, described in Eq. (3.2) is called the 
algebraic substitution technique, and converts rational 
polynomials in s into the desired rational polynomials in z. 

It will be shown that the selection of a specific f(z) 
to substitute for s, corresponds to the adoption of a numerical 
integration algorithm to integrate the phase state variables 
of the filter or system in discrete time. 



A. MEANING OF THE ALGEBRAIC SUBSTITUTION s = f(z) 

In order to find functions that can replace the variable 
s in the description of the continuous filter, H(s), the 
meaning of the substitution must be understood. In the first 
place, the variable s of the Laplace transform theory represents 
an operator in the time domain, i.e., it is a differentiator. 
Similarly, the inverse of s, s ^ , represents an integrator. 

Consider a filter whose transfer function is 

a + a. s + . . . + a s 

H(s) = — - , where m >_ n (3.3) 

b + b. s + . . . + b s m 
o 1 m 



If this transfer function is the ratio of the output signal 
transform to the input signal transform with zero initial 
conditions, then Eq. (3.3) can be expanded as follows 



= Y < s ) = W(s) . Y (s ) 
X (s) X (s) W(s) 



(3.4) 
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The following equations can be formed arbitrarily 



W(s) 



1 



(3.5) 



X(s) b Q + b^ s + b 



2 





m 




+ a s 
n 



n 



(3.6) 



Rearranging Eq. (3.5) and (3.6) yields 

m-1 

b s m W(s) = X (s) - E b. s 3 W(s) 
m j=0 3 



(3.7) 



H 

Y (s) = E a. s 3 W (s) 
j = 0 3 



(3.8) 



Equations (3.7) and (3.8) can be illustrated as in Fig. 

3.1, indicating a possible implementation, as is usually done 
in an analog computer simulation for low order systems. 

It has been shown that a transfer function of the type 
described in Eq. (3.3), can always be expanded in the form shown 
in Fig. 3.1. From this illustration it can be seen, that the 
replacement of the integrations (denoted by 1/s in the state 
space expansion of a continuous filter) by a function of z, 

1/f (z) , corresponds to a direct algebraic substitution o£ the 
variable s by the function s = f (z) in the closed form expression 
of Eq . (3.3). 

Furthermore, the linear operation l/f(z), should perform 
an approximate operation in the discrete time domain similar 
to the operation, 1/s, of the transform domain, if the responses 
obtained with the digital filter are to resemble the responses 
obtained from the continuous system. 
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Fig. 3.1, State Space expansion of a Transfer Function, 




In view of the above remarks, it can be concluded that 
the algebraic substitution technique represents the adoption 
of a numerical integration algorithm to perform the operation 
of integration for the state variables of the filter. 



B. OBTAINING RECURSIVE INTEGRATION ALGORITHMS 

To derive recursive integration algorithms, it must be 
decided whether or not the present value of the input is 
available to compute the present value of the output.. This 
question can be settled if some assumptions are made. For 
instance, if the time required to perform the computations is 
considered negligible, that is if the time elapsed between the 
present input and the corresponding output of the filter is of 
no importance, then the algorithms can be derived using inter- 
polation formulas. On the other hand, when this factor is 
considered important and the output at the present time is 
computed without the present input value, then some extrapola- 
tion formula must be used. This latter condition involves some 
form of a predicting method, and intuitively seems to be less 
accurate. Among several different criteria to find recursive 
integration algorithms, the following methods are described 
because they lead to well known integration formulas which are 
used in practice. 

a. The Adams-Moulton Method. 

This method is derived from the expression 



fti 



y k = y k-l 



f (t) dt 
x 



'k-1 



( 3 . 9 ) 
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“ s t 

where f (t) is a polynomial of (n-1) degree approximating 
the input signal, based on the last n available input values. 

The recursive integration algorithms for the first 
two degrees of approximating polynomials are 

Interpolating Form Extrapolating Form 

n (Using last input) (Not using last input) 

1 y k =y k-l f 2 (X k +X k-l ) y k =y k-l + J (3X k-l _X k-2 ) 



2 y k y k-l + 12 (5X k +8X k-l X k-2* y k y k-l + 12 (23X k-l 16X k-2 +5X k-3 * 



From these algorithms it can be seen that the first 
degree approximation using the last input value is the well 
known Trapezoidal rule of integration, 
b. The Milne's Method. 

This method correspond to the numerical integration 

criterion 



y k 



y k-n 




(t) dt 



(3.10) 



s t 

where f x (t) is the polynomial of (n-1) degree approximating 
the input function using the last n available values of the 
inputs . 



Clearly, when the approximating degree n is unity 
Adams-Moulton and Milne's methods are the same, giving identical 
recursive integration algorithms. For a polynomial approxima- 
tion of second degree the recursive Milne's algorithms are: 
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Interpolating Form 
n (Using last input) 



Extrapolating Form 
(Not using last input) 




+j (X k +4X 





It can be seen that Milne's first degree formula also leads to 
the trapezoidal rule and that the second degree formula leads 
to the Simpson's 1/3 rule of integration, when the last input 
is used. 

These recursive integration algorithms as well as others 
-have been analyzed in the literature under the general category 
of digital integrators, as for example in one of the first 
papers on the subject which was published by John M. Salzer 
in 1953 [3], or in a paper by A. P. Sage and R. W. Burt [11]. 

To obtain the frequency spectra of the algorithms, one 
must first obtain the equivalent Z transform representation. 

For example consider the trapezoidal rule of integration 



The equivalent digital filter can be derived by taking 
the Z transform 




(3. .11) 



Y (z) = z 1 Y ( z) + j (X (z) + z 1 X (z) ) 



( 3 . . 12 ) 



and the filter transfer function is 



T 1 + z 



-1 



-1 



(1.13) 



X (z) 



2 1 - z 
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This function corresponds to the algebraic substitution 



2 1 - z 

T 1 + z 



-1 



(3.14) 



The frequency spectrum of the filter described in (3.13) 
is obtained by substituting z = exp(jtoT) 



H (e > T ) =(£) L+_e 

2 1 



-jwT 



_ e - j wT 



(3.15) 



H (e 

and hence 



jwT) . t uoT 

= 2 c ° tan — 



13.16) 



|H(e^ wT ) | = cotan 
Arg H (e^ wT ) = - — 



(3.17) 

(3.18) 



These equations can be compared with the magnitude and 
phase of an ideal integrator 



H(s) 


1 

s 








(3.19) 


H(jco) 


1_ 

ju 








(3.20) 


H(jw)| = 


1 

w 








(3.21) 


H(jw) 


u 

2 








(3.22) 


Equation 


(3.17) 


can 


be expanded in a series 


of the 


form 


|H(e jwT ) | 


= 1 1 


( 2 _ 

k 0)T 


coT _ o) 3 T 3 _ \ ~ 1 _ 

6 360 J co 


coT 2 

12 


(3.23) 
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and therefore the magnitude spectrum of the trapezoidal rule 
digital integrator approximates the magnitude of an ideal 
integrator for values of co such that 



wT 2 

12 



<< 



1 

0 ) 



or equivalently 



0 ) << 



sn 

T 



or 



/3 ~ 1 

0) << — CO . - (0 

7T S 2 S 



(3.24) 



( 3 . 2 ! 5 ) 



(3.26) 



In Salzer's paper a study of the frequency spectra of 
several different digital integrators is considered with a 
comparison to the spectrum of an ideal integrator. 

Figure 3.2 shows curves of the trapezoidal rule and 
Simpson's 1/3 rd rule plotted against = co/co s . The straight 
line or real axis represents the ideal integrator and as can 
be seen the trapezoidal rule falls below this line, while 
Simpson's rule lies above this line. 

In general it can be said that the more complex the 
integration algorithm, the wider the band of coincidence with 
the ideal integrator and also the larger the number of multi- 
plications to be performed or the increase in computation 
time. As a result a possible figure of merit might be the ratio 
of 3 db bandwidth to the number of multiplications. As a matter 
of fact, simpler integration schemes perform as well or better 
than the complex ones in terms of the above figure of merit. 
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Fig. 3.2. Digital Integrators comparison, with 

ideal continuous Integrator, as function 
of normalized frequency. 
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C. THE BILINEAR SUBSTITUTION 

This particular case of the algebraic substitution 
technique corresponding to the adoption of the trapezoidal rule 
of integration has a different interpretation, as is shown by 
K. Steiglitz [5], and also discussed by Gold and Rader [15]. 

In order to avoid confusion between variables, the 
following notation will be used: 

s = 0 + jto (3.27) 



(6 + j-n.) T 
z = e J 



( 3 . 28 ) 



When the bilinear transformation is made, such that 



fl\ z - 1 
S T z + 1 



(3.29) 



then the entire jto axis is mapped onto the unit circle in the 
z plane. This can be shown from the inverse relationship of 



Eq. (3.29), namely 



T, 



z = 



+ <f> 



HD 



1 + tf ) s _ 

1 - (§) S 1 - (f) j<i> 



z = 1 



((2 ..30) 



To find the relationship between the variable to in the 
s plane and the variable in the z plane when s = jto,- the 
following equation may be used 

-1 cut 



(-) 

z = l ± = e j-^T = e j2 tan 

1 - (?) 



(3 .31) 



Equation (3.31) may be reduced to the form 



2 , toT 

■«■ = 5 ? arc tan — 



(3.32) 
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Equation (3.32) indicates the transformation relating the 
frequency axis of the s and z planes. Equation (3.29) means 
that when the bilinear transformation is applied to a rational 
expression in powers of s, H(s), a rational expression G(z) in 
powers of z is obtained. 



G(z) = H(| f-^-jj) (3.33) 

Furthermore, when the frequency spectrum of the s domain 
representation is obtained by evaluating the expression H(s) 
when s = j to , and H(jw) is obtained, and similarly if the fre- 
quency spectrum of the expression in powers of z is evaluated 
along the unit circle, the function obtained is 



G(e j " T) = H(§ 



arc tan 



out. 



(3.34) 



The resulting function has the same spectrum as H(jw) but 
compressed into the interval -tt/T < -n- <_ tt/T. 

The nonlinearity of the frequency scale can be overcome 
by prewarping [Ref. 15] or frequency scaling the chosen contin- 
uous realization, such that after the transformation the 
selected frequency is at the desired value. However, with this 
prewarping technique only one frequency may be scaled properly 
and thus certain characteristics of continuous filters as roll- 
off per octave are not preserved. This method of synthesis is 
particularly well suited to synthesize low pass filters with a 
cut off frequency much smaller than the sampling frequency, 
because then the transformation function for the frequencies, 
u) = — tan — ^ may be considered linear, preserving then thq 
characteristics of the continuous filter realization, from which the 
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characteristics were obtained. This situation can be forced 
by increasing the sampling frequency or equivalently by 
decreasing the sampling interval T, and this solution could be 
the cure for all problems, but this is not so. As was 
discussed in chapter two, reducing the sampling interval, is 
equivalent to reducing the time available for computations, 
and for a fixed computing speed, this results in a reduction 
in the maximum complexity of the filters that can be implemented. 
Therefore, the reduction in the sampling interval is the least 
desirable tool to be used. 

Besides this consideration of computing time, there is 
also another consideration about the number of bits required 
in the digital representation of the coefficients of the filter. 

First it must be noted that the coefficients in the digital 
filters have a finite accuracy, and hence the location of the 
poles in the z plane is made on a discrete basis. To clarify 
this consider a simple example where the location of a pole is 
directly related to the coefficient accuracy. The filter 

1 z" 1 

H (z) = — = - ^ (3.35) 

z - 1 - z 

has a pole in the z plane located at the point z = C^. Because 
of the finite accuracy in the representation of this pole 
can only be located at a finite number of positions. That is 
if a n bit binary representation of the coefficients is used, 
then the maximum accuracy in the pole position due to quantizing 
effect is 2 n . 
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Consider then the effect of a reduction in the sampling 
interval on the mapping of the s plane into the z plane under 
the bilinear transformation. 

To visualize this, a grid in the s plane, shown in Fig. 

3.3, is mapped into the z plane using the bilinear transforma- 
tion for values of T = 0.1, 0.05, 0.025, 0.125. The results 
are shown in Figs. 3.4, 3.5, 3.6 and 3.7. 

As can be seen, the area in the z plane into which, a 
given area of the s plane is mapped decreases with T, the 
sampling interval. Conversely a given area in the z plane 
represents a larger area in the s plane as T decreases. 

Consider a point z^ in the z plane distant one quantizing 
level from the point z = 1 (towards the origin on the real 
axis) and distant one quantizing level along a vertical direction 
(above the real axis) as shown in Fig. 3.8. 



= (1 “ ? n ) + j £ 



( 3 . 36 ) 



The corresponding point in the s plane is 




(3.37) 



S 1 T z.+l T 2-4 +j 4 
1 n J n 




(3.38) 




(3.39) 
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Fig. 3.3. Grid on the s plane to be mapped 

into the z plane using the bilinear 
transformation . 
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Fig. 3.4. Mapping of the grid of the s plane 
shown in Fig. 3.3, onto the z plane 
using the bilinear trnasformation 
for T = 0.1 sec. 
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Fig. 3.5. Mapping of the grid of the s plane 
shown in Fig. 3.3 onto the z plane 
using the bilinear transformation 
for T = 0.05 sec. 
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Fig. 3.6. Mapping of the grid of the s plane 
shown in Fig. 3.3 onto the z plane 
using the bilinear transformation 
for T = 0.025 sec. 
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Fig. 3.7. Mapping of the grid of the s plane 
shown in Fig. 3.3 onton the z plane 
using the bilinear transformation 
for T = 0.0125 sec. 
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Fig . 3.8. 



Point z = ( (l-€ n ) , jC n ) 
quantizing level from z 



distant one 
= 1 . 0 .. 
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2 g n (; n- 1)2 

T (2 - 5 n> 2 + 5 n 




(3.40) 



s 1 Z | + j where << 1 . 



(3 .41) 



The above derivation can be considered as the mapping cf 
a square of the z plane, each side measuring one quantizing 
level, £ , into an area in the s plane of side £ /T.. Because 
of quantizing effects a singularity at z^, can only be located 
at a finite number of positions, each of which corresponds to 
a position s^ . 

Equation (3.41) indicates that the grainniness of alloved 
locations in the s plane is dependent on the quantizing level 
and the sampling interval. When using a continuous filter 
design to synthesize a digital filter, the pole locations so 
obtained may not be conincident with these allowed positions, 
and so it may be necessary to shift them to those allowed 
locations which are closest to the design locations- Due to 
other considerations, like pole location sensitivity in the s 
plane, the maximum shifting allowable in the location of a pole 
in the s plane may be limited to a maximum distance A- If it 
is desired to make the maximum allowable shifting comparable 
to the grainniness in the s plane the following equation may 
be established 



A 



n_ 

T 



(3. .42a) 
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In the z -plane this grainniness would require n bits where 

2~ n = §n (3.42b) 

From (3.4 2) 

n = 3.3 log (^) (3.43) 

where n represents the number of bits used in the binary 
representation of coefficients and A is an indication of the 
grainniness of the grid in the s plane. 

The lower bound of the required number of bits can be 
determined from Eq. (3.43). Actually, the required number 
may be more, especially if the poles are located very far from 
the origin in the s plane, because in this case the distortion 
in mapping of areas is bigger and can no longer be neglected 
in computing A through this derivation. 

From the preceding discussion, it may be concluded that 
the finite representation of coefficients introduces an error 
in the location of the singularities of the function under- 
going the transformation, and that this error increases, that 
is represents a larger change in the s plane when the sampling 
interval is decreased. In other words, it is not possible to 
obtain a better approximation by decreasing only the sampling 
interval. Instead the number of bits which represent the 
accuracy of the coefficients must simultaneously be increased 
to obtain a better approximation to the analog specification. 

Because for real time operation all computations must be 
performed within a sampling interval, the above considerations 
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set a more rigorous limit than previously considered. This is 
so because an increase in bits to represent the coefficients 
tends to increase the computation time . 

D. CONCLUSIONS 

The algebraic substitution method provides an easy method 
to obtain a digital filter which approximates the character- 
istics of a continuous filter, providing then the desired link 
between the continuous filter theory and the discrete filter 
theory, with some restrictions. 

Of all the substitutions, the bilinear transformation is 

the easiest to apply and does not increase the complexity of 

th 

the filter, i.e., an order continuous filter is transformed 
th 

into a n*" order discrete filter. 

The simulation of continuous systems obtained by the 
substitution technique in general can be improved by decreasing 
the sampling interval. 

However, the spectrum obtained with a digital filter using 
this technique is only an approximation to the spectrum of the 
continuous filter from which it is derived, and if the con- 
tinuous filter is itself an approximation to a given specifica- 
tion there is no assurance that the digital filter will also 
be within the specifications. 
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IV. SYNTHESIS OF DIGITAL FILTERS 
FROM FREQUENCY SPECTRUM CHARACTERISTICS 



The standard Z transform technique and the algebraic 
substitution method considered in chapters two and three 
represent an attempt to extend the continuous filter theory to 

i 

discrete or digital filter theory. There is no method 
available which will guarantee that a digital filter will hcive 
an exact replica of the frequency spectrum of the continuous 
filter from which it is derived, except when the input signal 
and the filter transfer function are band limited and sampled 
according to Shannon's theorem using the Z transform technique. 
When using the above methods to approximate the frequency 
characteristics of the continuous filters (which themselves 
are an approximation to certain given specifications) there is 
no assurance that the resulting digital filter will meet design 
specifications also. 

A more direct approach to the study of the frequency 
spectra of digital filters seems to be necessary in order to 
provide better answers than the techniques derived from 
continuous filter theory. 

A part of the material presented here is an extension 

of an idea mentioned by M. Martin [6], in which a nonrecursive 

filter is constructed. The approximation of Martin to a given 

th 

specification using a Fourier series truncated at the n 
harmonic, requires 2 n words of storage and 2 n multiplications. 
The contribution described in this chapter leads to a recursive 
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filter, constructed from a ratio of two finite Fourier series 
approximating a given specification. This new procedure 
requires the same storage and number of multiplications, but 
allows one to obtain phase versus frequency characteristics 
other than the linear phase versus frequency with slope obtained 
when using Martin's method. 

The last part of the chapter presents the development of a 

generalized technique to obtain recursive digital filters 

derived from a magnitude squared criterion. The magnitude 

squared function is expanded as a ratio of two finite Rmrier 

th 

series. For the approximating series to involve up to the n 
harmonic, this procedure requires only n words of storage, 
which represents a substantial improvement over the magnitude 
methods discussed in the early part of the chapter. 

A. SYNTHESIS OF NONRECURSIVE FILTERS FROM A GIVEN MAGNITUDE 

SPECTRUM 

In certain applications only the magnitude spectrum of the 
output signal is important and a time delay of some sampling 
intervals can be tolerated. The procedure that follows allows 
the given magnitude function to be approximated with a finite 
summation of terms of a Fourier series. This procedure is 
mentioned by M. Martin [6] with reference to data transmission 
filters . 

Consider the elementary filter whose transfer function is 
of the form 

H (z) = | tl+z" 2 ] (4.1) 

A flow diagram of the filter is shown in Fig. 4.1a.. 
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(a) The elementary filter of first order 
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(b) The elementary filter of n th order 



Fig . 4.1. 



Elementary filters used for the 
Fourier series approximating technique 




The elementary filter described has a spectrum of the 



form 



tI , jwT. l r , -j2u)T. 

H(e J ) = j[l + e J ] 


(4.2) 


where z = exp(jcoT) has been substituted in Eq . (4.1) .. 

(4.2) can be reduced to the form 


Equation 


H(e^ U)T ) = cos uT e - ^' 1 


(4.3) 


and hence |H(e- ,uT ) | = |cos wT| 


(4.4) 


Arg H(e^ a> ^) = - coT 


(4.5) 



The elementary filter described in (4.1) is considered 
to be of order one and thus a family of elementary filters of 
higher order can be generated in the form of 



H n ( Z ) - | [1 + z- 2n ] 

Such filters have frequency spectra of the form 


(4.6) 


H n (e^ U)T ) = cos nwT e ^ ncoT 


(4.7) 


|H(e^ U)T ) | = |cos na)T| 


(4.8) 


Arg H(e- ,ujT ) = - nwT 


(4.9) 



The filter of order n is illustrated in Fig. 4. .lb.. 

Using these elementary filters it is possible to approxi- 
mate any given magnitude versus frequency characteristic 
|G(e- ,uT ) |, as will now be shown. 

First it must be noted that if the output of any of these 
filters described in Eq. (4.6) is delayed m sampling intervals, 
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the frequency response becomes of the form 




(4.10) 




cos ntoT e 



- j (n+m) toT 



(4.11) 



Comparing Eq . (4.11) with Eq . (4.7), it is clear that the 

magnitude of the frequency spectrum does not change by delaying 
the output. Only its phase changes in the form of an increase 
in the slope of the linear phase versus frequency 
characteristic . 

Therefore, by delaying the output signal of these filters 
a convenient number of sampling intervals, it is possible to 
make their output signals have the same phase characteristics, 
and hence if several of these specially delayed filters are 
added in parallel, the resulting magnitude spectrum is the sum 
of the individual magnitude spectra because all the complex 
quantities to be added have the same phase. 

An example of the above paragraph is illustrated in Fig. 
4.2a where the elementary filters of zero and first order are 
delayed such that their phase versus frequency characteristic 
become equal to the second order elementary filter. In Fig. 
4.2b the canonical form of the filter is shown, indicating 
symmetry in the coefficients. 

Generalizing, a filter synthesized as a summation of the 
delayed elementary filters in the form 
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(a) Equation form 




(b) Canonical form 



Y (z) 
X (z) 





+ 




z 



-4 



Fig. 4.2. Synthesis of nonrecursive filter 

approximating a given spectrum function 
involving up to the second harmonic of 
the Fourier Series. 
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(4.12) 



G (z) = 



N 
= Z 
n=0 



_ r , , -2n. . -N+n 

C [1 + z ] z 
n 



has a frequency spectrum which can be obtained by substituting 
= exp(ju)T) into Eq . (4.12), yielding 



G(e^ wT ) = E C cos na)T e ^ Na)T 
n=0 n 



(4.13) 



and hence |G(e^ wT ) | = 



N 

I E C cos nu)T 
1 n 

n=0 



(4.14) 



Arg G (e- ]ajT ) = - NojT 



(4.15) 



for real coefficients, C , with (n=0 , 1 , . . . ,N) . 

Equation (4.14) indicates the possibility of approximating 
any given magnitude spectrum with a Fourier series of finite 
numbers of terms. It is known that the magnitude spectra are 
even functions, and repetitive with period w s = 2 tt/T. This 
characteristic assures that it will always be possible to 
expand a magnitude spectrum function as a Fourier series with 
cosine terms only, assuring also the generality of Eq. (4.14). 

Equation (4.15) indicates that the resulting function 
will be delayed N sampling intervals, N being equal to the 
number of harmonics considered in the Fourier series expansion. 



If 



F 



n=0 



C cos nwT > 0 for all w 

n — 



(4.16) 



then the phase given by Eq. (4.15) will not display any dis- 
continuity jumps and will remain linear. This can always 
be accomplished by proper adjustement of the coefficient Cq . 
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It is well known that the Gibb's phenomenon occurs near 
a discontinuity of a function. In other words no matter how 
many harmonics are considered in the Fourier series approxima- 
tion, the summation of terms does not converge to the value 
of the function near a discontinuity, but does converge to ctn 
overshoot and a undershoot of approximately 10% of the 
discontinuous jump. In order to make certain that (4.16) is 
satisfied, the value of C Q must be adjusted to include the 
undershoot of the Gibb's phenomenon. In the case of an ideal 
low pass filter the Gibb's phenomenon is about 10% of the 
discontinuity. If Cq is adjusted so that the Fourier series 
is always positive, this implies that the attenuation in the 
stop band is approximately 20 db. 

Some other methods may be used to obtain a greater 
attenuation. For instance the Lanczos and Fejer methods to 
attenuate or eliminate the overshoot and ripple of the Fourier 
series can be used as discussed by R. W. Hamming [12] . 

A different approach to get around this inconvenience is 
to work numerically, with the given magnitude spectrum as a 
set of data points, through which a continuous function is to 
be fitted. The lack of discontinuities assures that Gibb's 
phenomenon will not occur. 

A simple example of nonrecursive digital filter synthesis 
is now presented. 

Synthesize a nonrecursive digital low pass filter approx- 
imating the magnitude function shown in Fig. 4.3, where w c =tt/ 3T 
and the order of the filter is chosen such that the approximation 
involves the fourth harmonic in the Fourier series expansion. 
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Representing the magnitude function as a Fourier series 



yields 



|G(e ja)T ) | = + E 



n=l 



C cos na>T 
n 



4 f U c 4 1 f W c 

where C = — cos ncoT dco = — — — sin na)T 

n to J co nT I 

s ' o s 1 o 



— — sin nTto 
nTco c 



4co sin nTco 
c c 



CO 



nTco 



For this example co /co = 1/6 and co T = tt/3. 

OS o 



Therefore the coefficients C become 

n 



C = — sin n tt/3 
n 3 n tt/3 



Cq = 0.66666 C 1 = 0.55133 C 2 = 0.27567 C 3 = 0.00 C 4 = -0. 

Consider first the case where C is not adjusted so that 
(4.16) is not satisfied. 



| G (e^ a)T ) | = 0.33333 + 0.55133 cos coT + 0.27567 cos 2coT 



+ 0.0 cos 3coT - 0.13783 cos 4coT 
and from Eqs. (4.14) and (4.12) 

, n -4 . 0.55133 ri , -2, -3 . Q..27567 

G(z) = 0.33333 z + - [1 + z ]z + = 



r , . -4, -2 0.13783 ri , -2 1 

[1 + z ]z - = [1 + z ] 



13783 
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and rearranging terms the filter becomes 



G (z ) = - 0.068915 + 0.13783 z 2 + 0.27567 z 3 + 0.33333 z ' 

+ 0.27567 z -5 + 0.13783 z“ 6 - 0.068915 z" 8 

In Figs. 4.4 and 4.5 the magnitude and phase versus 
frequency characteristics are shown for the filter of the 
example, where the coefficients are obtained directly from the 
Fourier series expansion. The discontinuities in the phase 
function can be seen to correspond to every cross over point 
shown in the magnitude function by an increase in attenuation. 
The lobes in the attenuation band of the magnitude function 
are the ripples of the Fourier series but distorted by the 
logarithmic scale. If C Q is adjusted so that C Q = 0.46030 then 

G (z) = - 0.068915 + 0.13783 z“ 2 + 0.27567 z -3 + 0.46030 z“ 

+ 0.27567 z -5 + 0.13783 z~ 6 - 0.068915 z“ 8 

The magnitude and phase response of this filter are shown 
in Figs. 4.6 and 4.7. 

In general the foregoing procedure requires 2 n words 
of storage for n harmonics considered in the Fourier series 
approximation. Furthermore the output signal has a delay of 
n sampling intervals due to the remarkable characteristic of 
this type of filters of having a linear phase with a slope 
dependent only in the number of delays of the filter, Failure 
to observe the condition stated in Eq. (4.16) produces discon- 
tinuities in the phase characteristic associated with each 
crossing of the w axis by the resulting Fourier series.. 
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(db) 



10 . 




-40. 



-50. 



0.0 0.2 0.4 0.6 0.8 1.0 

Fig. 4.4 Magnitude vs. normalized frequency. Non- 
recursive filter with coefficients from Fourier 
series . 
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Fig. 4.5 Phase vs. normalized frequency. Non- 
recursive filter with coefficients obtained 
from Fourier series. 
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Fig. 4.6 Magnitude vs. normalized frequency. Non- 
recursive filter, coefficients obtained 
from modified Fourier series to obtain 
linear phase. 
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Fig. 4.7 Phase vs. normalized frequency. Nonrecursive 
filter, coefficients obtained from modified 
Fourier series to obtain linear phase. 
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B. SYNTHESIS OF A SPECIAL TYPE OF RECURSIVE FILTER FROM 

A GIVEN MAGNITUDE SPECTRUM FUNCTION 

The method described in the previous section, where non- 
recursive filters were obtained for matching a given magnitude 
function, has been generalized by the author and allows a 
recursive filter to be obtained. 

This special type of filter has a Z transform of the 
impulse response 



H (z) 



, -1 

a + a , z + , 
n n-1 



-n -2n+l , -2n 

.+ a»z + . . .+ a , z + a z 

0 n-1 n 



b + b , z + . 
m m-1 



, , -m , . , -2m+l , , -2m 

.+ b z +...+ b ,z + b z 

o m-1 m 



(4.17) 



A canonical form of this filter is shown in Fig.. 4.8. 

Note that the polynomials of the numerator and the 
denominator have an even number of delays for the input and 
output signals, and furthermore are mirror image polynomials 
because of the symmetry in the coefficients. Such a filter 
can be considered to be derived from the expression 



H (z) 



A (z) 
B (z) 



(4.18) 



where A(z) = 



n 



. -1 , . -n , -2n+l -2n 

+ a ,z +...+ a rt z +...+ a z + a z 

n-1 0 n-1 n 



and B(z) = 



(4.19) 

, , , -1 , , . -m . . . -2m-fl , — 2m 

b + b .z +...+ b.z +...+ b .z + b z 

m m-1 0 m-1 m 



(4.20) 
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Fig. 4.8. Special type filter with 2 n delays for input 
and output signals and symetric coefficients. 





The frequency spectrum of such a filter is given by 




(4.21) 



and from Eqs . (4.12) and (4.13) 



n 



a + E 



2a, cos kwT 



o 



k=l 



- j (n-m) wT 




(4.22) 



m 

b + £ 2b, cos ktoT 

o , . k 



Equation (4.22) shows the properties of this special type 
of filter being considered. It indicates that its frequency 
spectrum is in the form of a ratio of two finite Fourier series, 
and the phase is linear with respect to frequency. 

When a magnitude versus frequency characteristic is given 
to be approximated with a digital filter, if an approximation 
in the form of the ratio of two finite Fourier series can be 
obtained, then a method to synthesize the filter is available 
by using Eq. (4.22) as will now be explained. 

1 . The Expansion of a Magnitude Function as a Ratio 
of Two Finite Fourier Series 

Consider first the properties of the Chebyshev 
polynomials, as discussed for instance by R. W. Hamming [13]. 



of both the Fourier series and the orthogonal polynomials; they 
are, in fact, the Fourier functions cos n 0 in the disguise 
of a simple transformation of the variable 

6 = arc cos x " 



I! 



The Chebyshev polynomials have all the properties 
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If a variable x is defined such that 



x = cos wT (4.23) 

from where -1 <_ X <_ 1 
then, 

T^(x) = cos (n arc cos x) = cos (wTri) (4.24) 

In view of Eq. (4.24) an expansion in Chebyshev 
polynomials in the variable x corresponds to a Fourier series 
in the variable ojT, when the two variables are related as in 
Eq. (4.23). This property can be used to approximate a given 
magnitude spectrum function as a ratio of two finite Fourier 
series in the following manner: 

Given a function |G(e-* wT ) | which has to be approxi- 
mated consider the following transformation 

|G(e^ (jT ) | = G-^ (u>T) = G^ (arc cos x) = H (x) (4.25) 

If this is done numerically, that is if pairs of 
values describing points of the magnitude versus frequency 
characteristic are given Eq. (4.25) can be considered as a 
change in the wT values of the given points, maintaining the 
ordinate values. The interest in doing this is that now the 
new function of x can be approximated as a ratio of two 
polynomials in x. Later the polynomials may be expressed as 
Chebyshev polynomials and then use of the property described 
in (4.24) leads to the Fourier series coefficients. That is 
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H(x) can be approximated by 



H (x) 



P n (x) 



(4.26) 



and numerator and denominator polynomials can be represented 
as an expansion of Chebyshev polynomials of the form 

n 

1 c k T k (x) 
k=0 K K 

H (x) : “ (4.27) 

E d T (x) 
k=0 K 



Using the property of Eq. (4.24) 
n 

E c. cos kwT 

k=0 k 

| G (e^^- 1 ) I : (4.28) 

m 

E d. cos kcoT 

k=0 k 

and Eq. (4.28) represents the approximation of a given function 
as the ratio of two finite Fourier series. 

For synthesis purposes, the coefficients in Eq. (4.28) 
can be equated with those of Eq. (4.22) yielding 

a k = c k k jt 0 (4.29) 

b k = d k k ^ 0 (4.30) 

The complete synthesis procedure is illustrated in 
Fig. 4.9. The following points should be noted: 
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Fig. 4.9. Procedure to find an approximating function 
to a given set of data points, in the form of 
the ratio of two finite Fourier series with 
cosine terms only. 
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a. If the discrete steps in the abcissae of 



|G(e^ wT ) | are uniformly spaced, the abcissae of H (x) will be 
nonuniformly spaced since x^ = cos T. 



b. Expressing H(x) = P (x) /Q (x) is possible in many 
ways. For the work presented here IBM subroutine ARAT is used. 

c. Expressing P (x) and Q(x) as Chebyshev poly- 
nomials is performed using the following tables taken from 
McCracken and Dorn [14]. 



1. Chebyshev Polynomials 




T 1 (x) = x 



T 2 (x) = 2 x 



2 



1 



T 3 (x) = 4 x 



3 



3 x 



(x) = 8 x^ - 78£ x 3 + 1 
T c (x) = 16 X 5 - 20 X 3 + 5 x 

D 



2. Powers of x as functions of T (x) 



n 



1 = T 



0 




TT(10T 1 +5T 3 +T 5 ) 



The degree of the numerator and denominator 



polynomials can be different, as indicated in Eq. (4.30) 



The 



difference in polynomial degree has direct influence on the 
phase characteristics, as shown in Eq . (4.22). Selecting equal 

degrees for numerator and denominator, it is possible to make 
the phase equal to zero for all frequencies. Otherwise the 
phase will change linearly with frequency, and the slope will 
be directly proportional to the difference in degree. 

An example of synthesis using this procedure is now 
presented. For comparison with the M. Martin method, the 
example given in the previous section is reworked. The 
magnitude spectrum of Fig. 4.3 will be approximated using a 
filter with eight delays, four in the numerator and four in the 
denominator, which implies the approximation of the given 
magnitude spectrum with a ratio of two Fourier series involving 
up to the second harmonic in each. The method used in the 
example to obtain the ratio of two Fourier series approximating 
the given magnitude spectrum function, is the numerical method 
illustrated in Fig. 4.9, using subroutine ARAT of the IBM 
scientific package to obtain the ratio of two polynomials in 
the transformed variable x = cos coT. The results are: 

Given eight points describing the magnitude spectrum 
function equally spaced on the wT axis in the range 

0 £ <JjT 1 11 (4 . 31) 

The change to the variable x is performed and a ratio of 
two polynomials in x approximating the points was found 

TT/ , 0.01184 + 0.02418 x + 0.01277 x 2 ,. A 

H (x) = ^ (4.32) 

1.000 - 1.39494 x + 0.47740 x 
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The polynomials are now expanded in terms of 
Chebyshev polynomials using table shown previously. This yields 



0.01823 T (x) + 0.02418 T (x) + 0.00639 T 0 (x) 

H (x) - ^ (4.33) 

1.23870 T q (x) - 1.39494 T ± (x) + 0.23870 T 2 (x) 

Using Eq . (4.24) the ratio of Fourier series can be 

obtained . 

| rf jwT. | _ 0.01823 + 0.02418 cos uT + 0.00639 cos 2wT 
1 ' ’ 1 1.23870 - 1.39494 cos U)T + 0.23870 cos 2wT 

(4 .34) 

From Eqs. (4.22) and (4.28) it follows that 

n , , 0. 00319+0. 01209z _1 +0.01823z _2 +0.01209z'' 3 +0.00319z“ 4 

G ( Z ) = 

0. 11935-0. 69747z +1.23870z -0.69747z +0.11935z 

(4.35) 

V 

The magnitude and phase of the spectrum has been 
plotted and illustrated in Figs. 4.10a and 4.10b. As a comparison 
with M. Martin's method shown in Fig. 4.4, it can be noted that 
there is an increase in attenuation and elimination of lobes 
as well as zero phase shift. 

C. SYNTHESIS OF RECURVSIVE FILTERS FROM MAGNTIDUE SQUARED 
FREQUENCY SPECTRUM FUNCTIONS 

In the previous section the filter is restricted to a 
minor image coefficient realization starting with the magnitude 
of a transfer function. In this section a new approach is 
developed to synthesize recursive digital filters, from a given 
magnitude squared function. 
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0.0 0.2 0.4 0.6 0.8 1. 



Fig. 4.10a Magnitude vs. normalized frequency. 

Recursive filter, coefficients obtained 
from numerical approximation described in 
Fig . 4.8. 
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0.0 0.2 0.4 0.6 0.8 1.0 

Fig. 4.10b Phase vs. normalized frequency. Recursive 

filter, coefficients obtained from numerical 
approximation described in Fig. 4.8. 
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It is shown that a digital filter, described in its most 
general form has a magnitude squared frequency spectrum which 
is a ratio of two finite Fourier series with cosine terms only. 
A method for obtaining an approximation to such a function, 
has been described in the previous section, using a change of 
variable and expansion of the approximating function in a 
series of Chebyshev polynomials. In this section the coeffi- 
cients of the filter are obtained by equating the coefficients 
of the Fourier series expansion with the coefficients of the 
filter according to some nonlinear relationships which are 
developed . 



1. 



filter , 
H (z) = 



Developement of the Theory 

Consider a generalized description of a digital 



of the form 

a o + a i Z_1 



+ a 2 z 



k 

-2 , -k E a.z 

+. . .+ a k z j = o 3 



“1 



m 



(4.36) 



—1 —2 — n 

b_ + b. z + b_ z +. . .+ b z E b. z 
0 12 m i=0 k 



-l 



The spectrum of this filter can be obtained by 
evaluating the magnitude and phase of Eq . (4.17) for values of 

z along the unit circle 

z = e^ wT = cos coT + j sin wT (4.37) 

Substituting (4.37) into (4.36), yields 



. a_+a 1 e^ a)T +a :) e^ ajT +...+ a v e 

H(a 3 “ T ) = -° .1 2 * 



b 0 + bl e-3“ T + b 2 e-J“ T 



+ . . . + b e 
m 



- jkwT 

- jmooT 



(4.38) 
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Consider 



( jwT = N_ ( j o)) 
} D ( jw) 



| H (e3“ T ) | 2 = 



D ( joo) 



H ( e jwT) | 2 _ Re 2 {N(jo))} + Im 2 {N(juj)} 
Re 2 {D(jw)} + Im 2 {D(jo))} 



(4.39) 



(a + Z cos nwT) + (Z a sin nu>T) 
|H(e^ T )| 2 = -° h=1 n=1 n 



m y m y 

(b + Z cos nwT) + (Z b sin nwT) 
o -i i -i n 

h=l n=l 



(4.40) 



which can be reduced to the form 



I jwT |2 = n=0_ n^O 

1 1 m „ m-1 



k 2 k-1 k-2 

Z a + 2Z a a , n cos ooT+2Z a a , ~ cos 2ooT +.. 
n „ n n+1 „ n n+2 

n=0 



m-2 

Z b ^+2Z b b cos ooT+2Z b b cos 2ooT +.. 

n n ~ n n+1 n n n+2 

n=0 n=0 n=0 



(4.41) 



|H( e j W T)|2 = 



Z c cos nooT 
n 



m 

Z d cos nu)T 

r» n 

n=0 



(4.42) 



Proof of Eq . (4.41) is given in Appendix A. Equation 

(4.41) indicates that the magnitude squared frequency spectrum 

is the ratio of two finite Fourier series with cosine terms 

only. Equation (4.42) shows the relationships between the 

coefficients of the Kmrier series, c and d and the coeffi- 

n n 

cients of the filter a^ and b^. Thus, from (4.22) and (4.23) 
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(4.43a) 



k-n 



n 



= 2 E 



j=0 



a . 
3 



j+n 



d 

n 



m-n 
2 E 
3=0 



b . 
3 




(4.43b) 



This relationship can be written in matrix form as 

follows 



a ^ 


a .. 


a _ . . . 


a. 








r* 


0 


1 


2 


k 




0 




c o 


0 


a^ 


a. . . . 


... a. 




a 


= 




0 


1 


k-1 




a i 




C l/2 


0 


0 


a„ . . . 


... a. 




a ~ 




c 






0 


k-2 




2 




2/2 


0 


0 


0 . . . 


. . . a„ 




a. 




c, , _ 








0 




k j 




_ c k/2 



Similar relationships hold for the denominator 

coefficients. As discussed in the previous section, it is 

possible to find an approximation to an even function as a 

rational expression of Fourier series so that c and d are 

n n 

known. The coefficients of the Fourier series, and hence the 
vector of c's in Eq. (4.44), will be obtained from the expansion 
of the desired spectrum function. Solving for the a's in 
Eq. (4.44) yields the coefficients of the numerator of the 
digital filter realization. An identical procedure can be 
followed for the coefficients of the denominator. 

This set of nonlinear equations has several solutions, 
in fact, as k gets larger, the number of solutions also 
increases. This is not surprising because of the nature of the 
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problem. That is, given a filter, a unique magnitude squared 
frequency spectrum function can be found, but the inverse 
procedure is not unique. Given the magnitude of the spectrum, 
the description of the filter is incomplete because there is 
no information about the phase characteristics. 

Consider an example of the use of this synthesis 
method to approximate a given magnitude squared function. For 
comparison with the example given the same magnitude spectrum 
shown in Fig. 4.3 will be squared. What results in the same 
function, which is now approximated with a filter using only 
two delays for the output and input signals. The nonlinear 
equations to be solved are of the form 



a 



0 



2 



+ 





(4.45) 



a l + a l 



a 2 



' 1/2 



(4.46) 



C 2/2 



(4.47) 



Solving this equations algebraically yields 




('4.48) 



a 



0 




a 



2 




(4.49) 



(4.50) 
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Equation (4.48) indicates there are four solutions 
for a^ and Eq. (4.49) indicates there are two solutions of a^ 
for each a^ . Hence there are eight sets of coefficients 
satisfying the given nonlinear equations. The same is true 
for the denominator set of coefficients and therefore there 
are 64 solutions for the filter. In general not all the 
solution may be valid if only real coefficients are being 
considered . 

Using the results obtained in previous example, 

Eq. (4.34), the magnitude squared function is expanded in the 
same ratio of Fourier series, yielding 

I , jwT, i2 _ 0.01823 + 0.02418 cos toT + 0.00639 cos 2wT ,-n 
“ 1.23870 - 1.39494 cos uiT + 0.23870 cos 2wT 1 ^ 

Substituting these coefficients into (4.48), (4.49), 

(4.50), two sets of solutions are obtained, called filter A 
and B respectively. 





FILTER A 


FILTER B 


a o 


0 . 038S83093 


0.081958605 


a l = 


0.099965521 


0.099965521 


a 2 


0.081958605 


0.038983093 


II 

o 

Xl 


-0.284898296 


-0.418921424 


b i = 


0.990978207 


0.990978207 


b, = 


-0.418921424 


-0.284898296 



The magnitude and phase versus frequency characteristics 
are shown in Fig. 4.11, 4.12, 4.13, 4.14. 
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Fig. 4.11 Magnitude vs. normalized frequency. 

Recursive filter, coefficients obtained 
from magnitdue squared criterion. Filter 
A. 
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Fig. 4.12 Phase vs. normalized frequency. Recursive 
filter, coefficients obtained from 
magnitude squared criterion. Filter A. 
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0.0 0.2 0.4 0.6 0.8 1.0 

Fig. 4.13 Magnitude vs. normalized frequency. 

Recursive filter, coefficients obtained 
from magnitude squared criterion. Filter 

B. 
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0.0 0.2 0.4 0.6 0.8 1.0 



Fig. 4.14 Phase vs. normalized frequency. Recursive 
filter, coefficients obtained from 
magnitude squared criterion. Filter B. 
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Comparison of filters A and B indicates the same 
magnitude for both filters and a different phase characteristic 
as was expected. With these sets of coefficients, still two 
more combinations are possible, that is the numerator of A 
and the denominator of B and vice versa. At this point, 
nothing has been said about the phase characteristic of this 
synthesis method. It is left as a suggestion for further 
research, the introduction of phase characteristics restrictions. 
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V. SUMMARY AND SUGGESTIONS 
FOR FURTHER RESEARCH 



The two main new ideas introduced in the thesis are: 

1. A generalization of Martin's synthesis technique to 
include recursive filters, which allows filters with zero 
phase or linear phase characteristics, to be obtained. 

2. A new synthesis technique, which depends upon the 
expansion of a given function as a ratio of two Fourier series. 
The resulting technique involves the solution of simultaneous 
nonlinear algebraic equations, but the filter requires half 
the storage of Martin's method for the same degree of approxi- 
mation. The number of multiplications to be performed is also 
decreased by one half. 

The foregoing studies have given rise to several basic 
questions for further research: 

1. Explore the concept of expanding a given function 
as a ratio of two finite Fourier series directly from Fourier 
series expansion theory, as versus the change of variable 
technique, involving Chebyshev polynomials, developed in this 
thesis. It would appear that an unlimited number of ratios 
are possible. If this is the case, filter stability can be 
assured by selecting the poles of the digital transfer function 
to lie inside the unit circle of the z plane, and then selecting 
the coefficients of the numerator to obtain the desired 
magnitude spectrum. 
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2. Consider synthesis procedures starting from a class- 
ical specification for the filter in terms of tolerance bands 
and attenuation rates. Determine how these design criteria 
can be met. 

3. Find a general method to solve the nonlinear equa- 
tions, Eq. (4.44), which were developed for the filter synthesis 
proceedure starting from a magnitude squared function as 
discussed . 
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Appendix A 

Derivation of Equation (4.22) 



k k 

|H(e- ,wT ) | 2 = (a n + E a cos ntoT) 2 +(E a sin nwT) 2 (A.l) 



0 , n 

n=l 



n=l 



n 



k j k 2 

b„ + E b cos ntoT)~+(E a sin mioT) 

0 n=l n n=l m 



consider the numerator only 



2 k „ k 2 

|N(ojT) | = a Q + E a cos nwT) z + (E a n sin nuT) 

n=l n n=l 

2 k k 2 ^ 2 

= a« +2a rt E a cos nwT+(E a cos nwT) + (E a sin nwT) 

0 0, n , n , n 

n=l n=l n=l 



(A. 2) 



The last two terms can be expanded as follows: 



k 2 k 2^2 2 

(E a cos nwT) + (E a sin nuT) = E a cos nooT 
, n , n , n 

n=l n=l n=l 



k-1 



2 . 2 



+ 2E a cos nwT E a. cos jwT + E a sin nooT 



n=l 



n 



j=n+l ^ 



n=l 



n 



k-1 



+ 2E a sin nwT E a. sin jcoT 
n=l n j=n+l - 1 



(A. 3) 



k-1 



^ 2 2 2 ^ x r 

E a (cos nwT+sin mooT)+2E a Icos nwT E a. 
n=l n n=l n 1 j=n+l 3 



cos jooT 



+ sin nwT E a. sin jwT;- 
n=n+l ^ > 



(A. 4) 
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2 



k-1 f k 



E a + 2 E <Z a cos nwT a. cos iuT 
n=l n n-1 X i=n+l n 3 



+ E a cos nuT a. sin juTl 



E a + 2 
*i n 
n=l 



n=n+l n 

k-1 k 

E 

j 



K ” _L K 

E { E a a. {cos nwT cos imT 
n-1 1 i=n+l n 3 



+ sin nwT sin jwT} 



I 



Use the trigonometric identity 

{cos nojT cos jwT + sin nuT sin jcoT} = cos (j-n)wT 



* 2 k-i k ^ 

Ea +2E xE aa. cos (j-n)ujT > 
i=l n n=l j=n+l n 3 3 



changing subscripts, i = j-n , j = n+i the limits of the 
second summation become 

j = n+1 i = 1 

j = k -> i = k-n 



n=l 



n 



k-1 


k-n 


E 


E a 


n=l 


i=l 


order Eq. 


k-1 


K-i 


E 


{ 2 


i=l 


L n=l 



n n+i 



Changing summation order Eq. (A. 8) becomes 
k 



a)T 



(A. 5) 



(A. 6) 



(A. 7) 



(A. 8) 



(A. 9) 
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Combining (A. 9) with (A. 2) yields 



2 2^ k 2 

|N(a)T) I = a n + 2 £ a.a., . cos i wT + E a 

0 i=l 0 0+1 n=l n 



k-1 k-1 

+ 2 £ E a a , . cos i wT 
. , , n n+i 

i=l n=l 



(A. 10) 



k 2 k-1 

= E a +2 a»a, cos k u)T + 2 E a^a„ , . cos i wT 
~ n Ok . . 0 O+i 

n=0 i=l 



k-1 k-1 

+ 2 E E a a . cos i u)T 
. , , n n+i 

l-l n=l 



(A. 11) 



k „ k-1 k-i 

= E a +2 a^a, cos k wT + 2 £ E a a . . cos i u)T 
„ n Ok . , „ n n+i 

n=0 i=l n=0 



k 2 k k-1 

= E a + 2 E E a a . . cos i wT 
» n . , ,, n n+i 

n=0 i=l n=0 



(A. 12) 
(A. 13) 



The same relationship holds for the denominator, and so 
established the validity of Eq. (4.22). 
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